Data Description:¶
The actual concrete compressive strength (MPa) for a given mixture under a specific age (days) was determined from laboratory. Data is in raw form (not scaled). The data has 8 quantitative input variables, and 1 quantitative output variable, and 1030 instances (observations).¶
Domain:¶
Cement manufacturing¶
Context:¶
Concrete is the most important material in civil engineering. The concrete compressive strength is a highly nonlinear function of age and ingredients. These ingredients include cement, blast furnace slag, fly ash, water, superplasticizer, coarse aggregate, and fine aggregate.¶
Attribute Information:¶
- Cement : measured in kg in a m3 mixture
- Blast : measured in kg in a m3 mixture
- Fly ash : measured in kg in a m3 mixture
- Water : measured in kg in a m3 mixture
- Superplasticizer : measured in kg in a m3 mixture
- Coarse Aggregate : measured in kg in a m3 mixture
- Fine Aggregate : measured in kg in a m3 mixture
- Age : day (1~365)
- Concrete compressive strength measured in MPa
Learning Outcomes:¶
- Exploratory Data Analysis
- Building ML models for regression
- Hyper parameter tuning
Import Libraries¶
In [423]:
# import libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
#%config InlineBackend.figure_format = 'retina'
from scipy.stats import zscore
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
# Import Linear Regression machine learning library
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
# Import KNN Regressor machine learning library
from sklearn.neighbors import KNeighborsRegressor
# Import Decision Tree Regressor machine learning library
from sklearn.tree import DecisionTreeRegressor
# Import ensemble machine learning library
from sklearn.ensemble import (RandomForestRegressor, GradientBoostingRegressor,AdaBoostRegressor,BaggingRegressor)
# Import support vector regressor machine learning library
from sklearn.svm import SVR
#Import the metrics
from sklearn import metrics
#Import the Voting regressor for Ensemble
from sklearn.ensemble import VotingRegressor
# Import stats from scipy
from scipy import stats
#importing the metrics
from sklearn.metrics import mean_squared_error,mean_absolute_error,r2_score
#importing the K fold
from sklearn.model_selection import KFold
#importing the cross validation score
from sklearn.model_selection import cross_val_score
#importing the preprocessing library
from sklearn import preprocessing
# importing the Polynomial features
from sklearn.preprocessing import PolynomialFeatures
#importing kmeans clustering library
from sklearn.cluster import KMeans
from sklearn.utils import resample
Import Data¶
In [424]:
# import data
df=pd.read_csv("concrete (1).csv")
df.shape
Out[424]:
Data Types and Data set values Analysis¶
In [425]:
# get info on data avialable and compare it with details provided about data set
df.info()
- Total 1030 rows and 9 columns present
- All data are numerical as expected from the description
- As all values are numeric then there is less chances of data having text or symbols apart from few unexpected numbers if any as outliers
- Non- NULL count is same for all column suggests no null values
- We have 8 independent and 1 dependent values
- Target Column: -Strength
In [426]:
#Analyze distribution
df.describe().T
Out[426]:
- All columns are qantitative
- There is huge gap between min and max values for columns
- Mean and 50% are not similar for many cloumns, represents skewness for those data.
- Mean is higher for slag, ash,age compared to 50%, represents right skewness (long tail towards right)
- many min values are at 0, need to validate its feasibility and practicality
In [427]:
#view the top 5 rows of data
df.head()
Out[427]:
In [428]:
#view bottom 5 rows of data
df.tail()
Out[428]:
In [429]:
# view random 10 data from dataset
df.sample(10)
Out[429]:
- With head, tail and random samples we can notice there is a diverse spread of data over all attributes.
Analyzing Attributes¶
In [430]:
fig, axs = plt.subplots(nrows = 9, ncols=2, figsize = (10,30))
for i, x in enumerate(df.columns):
sns.distplot(df[x],ax=axs[i,0],bins=50, rug=True)
sns.boxplot(orient='v', data=df[x], ax=axs[i, 1])
fig.tight_layout()
plt.show()
plt.clf()
plt.close()
#for colname in df.columns:
# sns.distplot(df[colname])
- Cement: has close to normal distribution (slight tail towards higher denominations of numbers) and no visual outliers present, mode seems to be around 200-250
- Slag: mode of the distribution is clearly 0, rest of the dataset seems like right skewed as it has few higer numbers beyond 300. Few outliers alos visible beyond 350.
- Ash: Similar to Slag, this also has a mode of 0, but inlike Slag it has a huge gap of values between 0 and 50 or 100. This suggests the value 0 might need a little more investigation to deicde if its a missing value or something else. No visible outliers noticed.
- Water: ditributions has multiple peaks but its look spreaded across both side of peak. Few outliers present in both lower and higher side of values
- Superplastic: mode is zero and its right skewed. few outlier present in higher side of values
- Coarseagg: datapoints looks distributed and no outlier presence
- Fineagg: Distribution seems normal but has few outlirers on higher values
- Age: highly right skwed data points, lot of outlier presence
- Strength: Data points seems distributed normally but few outlier presence towards higher values
In [ ]:
In [431]:
df.boxplot(figsize=(15, 10));
- From box plot we can clearly see there is high presence of outliers on slag, water, superplastic, fineagg, age and strength
In [432]:
sns.pairplot(df, diag_kind='kde');
- From intital view its difficult to see any strong liner or polynomial correlation among the attributes. We can observe ther are few week correlation present between few attributes. Crrelation heatmap shold provide more insights.
- Majority of bi-variate graphs shows a cloud like structure
- From diagonal analysis on kde plots, we can observe there is presence of minimum two distinct peaks in the independent attributes
In [433]:
df.hist(bins=30,figsize=(25, 10));
In [434]:
plt.figure(figsize=(15, 10))
sns.heatmap(df.corr(),annot=True,cmap="cividis",linecolor="black" );
In [435]:
df.corr()['strength'].sort_values()
Out[435]:
- From heat map we can clearly observe the linear relations.
- Majority of attributes has week correlations ranging from -3 to +3.
- Few strong -ve correlation are visible among attributes like between: superplastic & water, superplastic & ash, fineagg & water
- from correlation details with strength we can notice highest positive/negative correlation is with cement, superplastic, age and water.
- We can plan for removal of low correlated attributes down the line based on analysis.
In [436]:
sns.lmplot(x="cement",y="strength",data=df)
plt.show()
- we can notice a clear positive correlation
In [437]:
#cement vs water
sns.lmplot(x="cement",y="water",data=df)
plt.show()
- We cal clearly see the line is almost parallel to x axis, representing very minimal correlation between the attributes.
In [438]:
sns.lmplot(x="superplastic",y="water",data=df)
plt.show()
- we can notice the strong negative correlation line
In [ ]:
In [439]:
sns.lmplot(x="ash",y="strength",data=df)
plt.show()
- we can notice the line is almost paralle to the x axis.
Missing Values¶
In [440]:
df.isnull().sum()
Out[440]:
- No null values present
Analyzing Outliers and treatment¶
In [ ]:
In [441]:
df.boxplot(figsize=(10, 10));
- From box plot we can clearly see there is high presence of outliers on slag, water, superplastic, fineagg, age and strength
In [442]:
# NUmber of outliers
In [443]:
for cols in df.columns[:-1]:
q1=df[cols].quantile(0.25)
q3=df[cols].quantile(0.75)
iqr=q3-q1
low=q1-1.5*iqr
high=q3+1.5*iqr
print('Outliers count for',cols, df.loc[((df[cols]>high) | (df[cols]<low)),cols].count())
print('high value Outliers count for',cols, df.loc[(df[cols]>high),cols].count())
print('low value Outliers count for',cols, df.loc[(df[cols]<low),cols].count())
- there is significant number of outliers presnt in age attribute compared to other attributes.
In [444]:
#Records containing outliers
for cols in df.columns[:-1]:
q1=df[cols].quantile(0.25)
q3=df[cols].quantile(0.75)
iqr=q3-q1
low=q1-1.5*iqr
high=q3+1.5*iqr
print()
print('Outliers count for',cols)
print(df.loc[((df[cols]>high) | (df[cols]<low)),:])
- we will be treating the columns with same logic used in box plot using IQR
- We are replacing the Outliers with boundry values to keep the higer values colse to higher and lower boundries.(Capping)
- Goal is to keep the value inside acceptable ranges without loosing its main charatcteristic of being outlier or being a higher/lower number. Using this approach we will slightly increase the frequency on boundry values but they will be inside range. We have to keep an eye ot for a bigger group formation towards the boundry or a peak in graphs.
- Replacing with Mean or median might miss represent the main characteristics of the data set which is pushing it to be be outlier.(Introducing bias by making it a normal curve)
- We will be doing this only once, we may see new outliers as the IQR might change
In [445]:
#- selecting all but leaving out the last columns whihc is our target
for cols in df.columns[:-1]:
q1=df[cols].quantile(0.25)
q3=df[cols].quantile(0.75)
iqr=q3-q1
low=q1-1.5*iqr
high=q3+1.5*iqr
df.loc[(df[cols]<low),cols]=low
df.loc[(df[cols]>high),cols]=high
In [ ]:
In [446]:
#Rechecking the outliers and qartiles
df.boxplot(figsize=(10, 8));
- No more outliers present in the independent columns
In [447]:
#to capture if we are getting any additional peaks or bulges on the due to outlier treatments
sns.pairplot(df, diag_kind='kde');
- Peaks are almost similar to what we have noticed earlier which suggestes our outlier treatments has not affected highly towards the borders.
In [ ]:
In [ ]:
Feature Engineering, Model Building and Model Tuning¶
- All the feature showed less correlation internally,
- few features have very less correlation with Strength as well which we can identify and eliminate
- we will be using lasso to identify which all features have 0 contributing factor and
- compare with what we have analyzed during correlation heatmap analysis
In [449]:
#Scaling data set
df_scaled = preprocessing.scale(df)
df_scaled=pd.DataFrame(df_scaled,columns=df.columns)
In [450]:
df_scaled.describe().T
Out[450]:
In [451]:
#spliting dependent and independent variable
X=df_scaled.iloc[:,0:8]
y = df_scaled.iloc[:,8]
In [452]:
# Training and test set in 70:30 ratio
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.3, random_state = 50)
In [ ]:
In [ ]:
In [453]:
# OLS - LinearRegression
# We have noticed earlier that all the features were not very strong predictors of strength (max corr was 66% with cement)
# we are expecting our Linear model to be a weak predictor
regression_model = LinearRegression()
regression_model.fit(X_train, y_train)
for idx, col_name in enumerate(X_train.columns):
print("The coefficient for {} is {}".format(col_name, regression_model.coef_[idx]))
intercept = regression_model.intercept_
print("\nThe intercept for our model is {}".format(intercept))
In [454]:
# predict mileage (mpg) for a set of attributes not in the training or test set
y_pred = regression_model.predict(X_test)
# Since this is regression, plot the predicted y value vs actual y values for the test data
# A good model's prediction will be close to actual leading to high R and R2 values
plt.scatter(y_test, y_pred);
- We can notice the actual and predicted values are spread widely and not concentrated towards the diagonal
In [ ]:
In [455]:
#Ridge
ridge = Ridge(alpha=.3)
ridge.fit(X_train,y_train)
#print ("Ridge model:", (ridge.coef_))
for idx, col_name in enumerate(X_train.columns):
print("The coefficient for {} is {}".format(col_name, ridge.coef_[idx]))
In [456]:
#Lasso
lasso = Lasso(alpha=0.1)
lasso.fit(X_train,y_train)
#print ("Lasso model:", (lasso.coef_))
featureAnalysisLasso=lasso.coef_
for idx, col_name in enumerate(X_train.columns):
print("The coefficient for {} is {}".format(col_name, featureAnalysisLasso[idx]))
#we will use the lasso coeff details later as well for feature selecction
In [ ]:
In [457]:
#Comparing all linear regression details
In [458]:
print(regression_model.score(X_train, y_train))
print(regression_model.score(X_test, y_test))
In [459]:
print(ridge.score(X_train, y_train))
print(ridge.score(X_test, y_test))
In [460]:
print(lasso.score(X_train, y_train))
print(lasso.score(X_test, y_test))
In [ ]:
- As expected from linear correlation details, OLS(linear regression) is not a strong predictor and able provide 71.5% accuracy with Train and 75.2% with test set.
- With Ridge we can notice the results are very close to OLS
- With Lasso we can notice a significant drop in accuracy as it has dropped 3 features from its model consideration (fineagg, coarseagg, ash)
- we will try with polynomial approach to see if there is a significant cahnge in it.
- From the graph analysis we haven't noticed any variable which shows a polynomial relation as well but we will git it a try.
In [ ]:
In [ ]:
In [461]:
#Polynomial approach and complexity
In [462]:
from sklearn.preprocessing import PolynomialFeatures
#going with degree 2 only
poly = PolynomialFeatures(degree = 2, interaction_only=True)
#poly = PolynomialFeatures(2)
#X is already scaled and will use it
X_poly = poly.fit_transform(X)
X_train, X_test, y_train, y_test = train_test_split(X_poly, y, test_size=0.30, random_state=50)
X_train.shape
Out[462]:
- Shape of train data set changed to 721, 37
- we have 37 cloumns now instead of 8 (close to 5 fold increase in variables)
In [463]:
regression_model.fit(X_train, y_train)
print(regression_model.coef_)
In [464]:
print(regression_model.score(X_train, y_train))
print(regression_model.score(X_test, y_test))
In [465]:
ridge = Ridge(alpha=.3)
ridge.fit(X_train,y_train)
print ("Ridge model:", (ridge.coef_))
In [466]:
print(ridge.score(X_train, y_train))
print(ridge.score(X_test, y_test))
In [467]:
lasso = Lasso(alpha=0.01)
lasso.fit(X_train,y_train)
print ("Lasso model:", (lasso.coef_))
In [468]:
print(lasso.score(X_train, y_train))
print(lasso.score(X_test, y_test))
- Comparing the results of polynomial features with linear, ridge and Lasso, we can notice there is significant increase in the accuracy for the traina nd test data set.
- For Lasso we can notice a significant jump in accuracy from 63% to 78% (eliminating 12 out of 37 features).
- Polynomial features are able to dig out relations which were not visible in the data analysis phase but we also observed multi fold increase in features which also increases the complexity of the model. We may see more improvement by increaseing the degree of the polynomial features but the complexity will be too high then. We are ending our analysis with polynomial features here and will try out few other models to see if we are able to get better numbers.
In [ ]:
Explore for gaussians. If data is likely to be a mix of gaussians, explore individual clusters and present your findings in terms of the independent attributes and their suitability to predict strength¶
In [469]:
sns.pairplot(df_scaled,diag_kind='kde')
Out[469]:
- We can clearly notice there are minimum two gausian or clear peaks visible in the data sets (slag, ash, superplastic, age). WE max observe max of 3 peaks for Age but its comapratively smaller then others.
- Only field which has clear linear relationship with strength is cement, for all other feture we can see the distibution is spread across strength on y axis for any given point of feature value in x.
- even in the cloud formation in the graphs, we can observe there is two clear distinct groups froming one which is close to zero and rest others
- Based on above pairplot and below correlation details with strength (dependent) column
- Cement: -- Cement seems to have clear and maximum positive correlation, data distribution looks close to normal. This feature should be included in the model preparation.
- Age: -- For age we can clearly notice formation of groups for strength at diff age levels and is a positive correlation. Seems to have multipe gaussians present.This feature should be included in the model preparation.
- superplastic: -- superplastic seems to have positive correlation but comaparatively less than cement or age from graphs, seems to have two gaussian or groups or peaks present from the graphs analysis. This feature should be included in the model preparation.
- slag: -- slag seem to have a weak relation with strenght. There is clear indication of two gaussian. While comparing with strength we can notice there is two group formation one less than 0 and other above it. This has some significance but this feature is on the list for elimination. We will finalize the details based on further analysis.
- ash: -- This feature has clearly weak correlation with strength. Two gaussinas are clealry visible. On the scatter plot with strength we can clearly see two distinct cloud formations but both doesnt have a clear head or tail. This is feature will not be beneficial in the model preparation.
- coarseagg,fineagg: -- both these features have similar data distribution and correlation values. From Correlation details with strength we can notice they have negative correlation and are better than ash but visually its very hard to identify any head or tail for it. This is feature will not be much beneficial in the model preparation but will keep it for further analysis as they have higher correlation value than ash and slag.
- water: -- water has high negative correlation comapred to other features having negative correlation. We can notice 3 gaussians in the distribution. On the scatter plot with strength we can notice thre group formation on less than -2, more than +2 and in between -2 to +2. This feature should be included in the model preparation.
In [470]:
df_scaled.corr().strength.sort_values()
Out[470]:
In [ ]:
In [471]:
# Trying to see how many groups we can capture with kmeans & GaussianMixture models
from sklearn.cluster import KMeans
cluster_range=range(1,15)
cluster_error=[]
for num_clusters in cluster_range:
cluster=KMeans(num_clusters,n_init=20)
cluster.fit(temp)
labels=cluster.labels_
centroid=cluster.cluster_centers_
cluster_error.append(cluster.inertia_)
clusters_df = pd.DataFrame({"num_clusters": cluster_range, "cluster_errors": cluster_error})
clusters_df[0:15]
from matplotlib import cm
import matplotlib.pyplot as plt
plt.figure(figsize=(6,4))
plt.plot(clusters_df.num_clusters, clusters_df.cluster_errors, marker = "X")
Out[471]:
- Kmeans shows a weak prediction (a curve), which suggests its not able to identify groups in data sets clearly
- but from kmeans clustering we can see the elbow formation at 2
In [ ]:
In [ ]:
In [472]:
from copy import deepcopy
temp=deepcopy(df_scaled.drop(columns=['strength'], axis=1))
temp2=deepcopy(df_scaled)
In [473]:
#Splitting the data in 2 set
#pandas will cut the data in 2 ranges and convert the continuous data to categorical for this purpose
val2=pd.cut(temp2.strength, bins=2, labels=np.arange(2), right=False)
In [474]:
# training gaussian mixture model
from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(n_components=2, random_state=50)
gmm.fit(temp)
Out[474]:
In [475]:
#predictions from gmm
labels = gmm.predict(temp)
frame = pd.DataFrame(temp)
frame['cluster'] = labels
frame['stregth_cluster']=val2
frame.columns = ['cement', 'slag', 'ash', 'water', 'superplastic', 'coarseagg',
'fineagg', 'age', 'cluster','stregth_cluster']
color=['blue','green','cyan', 'black','red']
for k in range(0,2):
data = frame[frame["cluster"]==k]
plt.scatter(data["cement"],data["water"],c=color[k])
plt.show()
#plotting the clusters to identify if cluters created from modela nd clusters created in data are any match
#plot with split as per model prediction
In [476]:
for k in range(0,2):
data = frame[frame["stregth_cluster"]==k]
plt.scatter(data["cement"],data["water"],c=color[k])
plt.show()
#Plot with real split with strength
- we can see there is a significant similarity on identifying the dots towards far right by model but with the cloud we can see a mix of all data.
In [477]:
frame.head(10)
Out[477]:
In [478]:
from IPython.display import display
In [479]:
ct = pd.crosstab(frame['cluster'], frame['stregth_cluster'])
display(ct)
- From cross tab we can see there is a clear mismatch in groups and there are significant laps in recall and preceission.
In [ ]:
In [ ]:
Feature Importance and feature selection¶
In [480]:
#Feature selection
from mlxtend.feature_selection import SequentialFeatureSelector as sfs
In [481]:
# Build Lin Reg to use in feature selection
linR = LinearRegression()
# Training and test set in 70:30 ratio
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.3, random_state = 50)
# Build step forward feature selection
# we have choosen 5 as from our earlier exploration with lasso and feature analysis
sfs1 = sfs(linR, k_features=5, forward=True, scoring='r2', cv=5)
# Perform SFFS
sfs1 = sfs1.fit(X_train.values, y_train.values)
In [482]:
sfs1.get_metric_dict()
Out[482]:
In [483]:
from mlxtend.plotting import plot_sequential_feature_selection as plot_sfs
import matplotlib.pyplot as plt
fig = plot_sfs(sfs1.get_metric_dict())
plt.title('Sequential Forward Selection (w. R^2)')
plt.grid()
plt.show()
In [484]:
# Which features?
columnList = list(X_train.columns)
feat_cols = list(sfs1.k_feature_idx_)
print(feat_cols)
In [485]:
subsetColumnList = [columnList[i] for i in feat_cols]
print(subsetColumnList)
- Above is the list of features for model selection based on importance.
In [486]:
#Lasso details for feature
for idx, col_name in enumerate(X_train.columns):
print("The coefficient for {} is {}".format(col_name, featureAnalysisLasso[idx]))
- From details from Lasso feature selection, SequentialFeatureSelector and correlation Heat map we can clearly identify the 5 features we can go with for further model building and testing without loosing too much data.
- Features identifed for Model - 'cement', 'slag', 'water', 'superplastic', 'age'
In [ ]:
In [487]:
# Split X and y into training and test set in 70:30 ratio
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.3, random_state = 50)
dt_model = DecisionTreeRegressor()
dt_model.fit(X_train , y_train)
Out[487]:
In [488]:
#printing the feature importance
print('Feature importances: \n',pd.DataFrame(dt_model.feature_importances_,columns=['Imp'],index=X_train.columns))
- From decission tree regressor as well we can notice we are ending with same features, which are cement,slag,water,superplastic, age
Model Building¶
- During Linear correlation analysis, we have noticed there is very less attributes in data set which have high correlation or pridicting power for target attribute.
- During our linear and polynomial analysis we have noticed the accuracy are not too high.
- Polynomial feature were able to provide slightly improved results compared to linear but with high number of attributes.
- Support vector regressor or KNN are strong models for classification, they have significant capabilities for regression problems as welll but doesnt seems like a suitable candiate for this regression predictions.
- Considering above details, we will put our model building efforts on ensamble methods
- We are proceeding with those attributes only which have significant contribution on targets predictions
In [578]:
modelComp=pd.DataFrame()
In [570]:
X=df_scaled.iloc[:,0:8]
y = df_scaled.iloc[:,8]
seed=50
num_folds = 50
#Removing less contributing features
X=X.drop(['ash','coarseagg','fineagg'],axis=1)
# Split X and y into training and test set in 70:30 ratio
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.3, random_state = seed)
In [491]:
X.columns
Out[491]:
In [492]:
#Decission tree Regreesor with default parameters and defaul values
dt_model = DecisionTreeRegressor()
dt_model.fit(X_train , y_train)
y_pred = dt_model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',dt_model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',dt_model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
In [493]:
from scipy.stats import pearsonr
sns.jointplot(x=y_test, y=y_pred, stat_func=pearsonr,kind="reg");
- For decission tree with default sets and we can see that training set accuracy is close 99.5% but test set accuracy is 84.7%. This variation suggests we have an overfitting model.
- with pearsonr plotting as well we can see the values are spread far from the expected line.
In [494]:
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
model1 = dt_model
results = cross_val_score(model1, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
- This model with default set is not a suitable model as it has 77% average accuracy but having a 23% variance, from detils we can notice it predicts values as low as 60% and for few data sets its around 95%. This suggests the model is overfitted and will trade poorely with many data sets.
Prunning Decission tree¶
In [579]:
#Regularizing Decission Tree with diff values
model = DecisionTreeRegressor( max_depth = 8,random_state=seed,min_samples_leaf=4)
model.fit(X_train , y_train)
Out[579]:
In [580]:
y_pred = model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
#model1 = model
results = cross_val_score(model, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
modelComp=modelComp.append(pd.DataFrame({'Model':['Decission Tree'],
'Train Accuracy':[model.score(X_train,y_train)],
'Test Accuracy':[model.score(X_test , y_test)],
'Kfold-Mean-Accuracy':[results.mean()],
'Kfold-StdDeviation':[results.std()]}))
- With Prunned decission tree we can see the Train and test data sets have close by accuracies. This reduces our Overfitting
- But considering the cross val score and details - we can notice it still has very high std. deviation. But there is still a high gap between the train and test accuracies and we are still ending up with overfitting model. Any further prunning of the tables leads to reduced test and cross val accuracies.
- We will move ahead with some other models and will try for Grid search on the models which have better performances in next steps.
Random Forest¶
In [497]:
#Random FOrest With default config
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor()
model.fit(X_train, y_train)
Out[497]:
In [498]:
y_pred = model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
#model1 = model
results = cross_val_score(model, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
- With default values only we are able to get a significant better value then Decission tree
- Train Data set is close to 98% accuracy which hints towards overfit but test data set is as well around 89.8%, which is not very far from train accuracy
- Model seems slightly over fit but we can give it a try with grid search to see if its make the result any better and try to close the gap between train and test accuracies.
- Grid search on Random Forest as its default is better than Decission tree
In [499]:
from sklearn.model_selection import GridSearchCV
parameters = {'bootstrap': [True],
'max_depth': [5, 10, 15, 20, 25, 30, 35, 40, 50],
'max_features': ['auto', 'sqrt'],
'min_samples_leaf': [1, 2, 4, 8],
'n_estimators': [100]}
clf = GridSearchCV(RandomForestRegressor(),
parameters,
cv = 5,
verbose = 2,
n_jobs= 4)
clf.fit(X_train, y_train)
clf.best_params_
Out[499]:
In [582]:
#Using the grid search is giving us better train data accuracy but its leading it to overfit zone
#so prunned the max depth after few iterations
model = RandomForestRegressor(bootstrap = True,
max_depth= 9,
max_features= 'sqrt',
min_samples_leaf= 1,
n_estimators= 100)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
#model1 = model
results = cross_val_score(model, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
modelComp=modelComp.append(pd.DataFrame({'Model':['Random Forest Regressor'],
'Train Accuracy':[model.score(X_train,y_train)],
'Test Accuracy':[model.score(X_test , y_test)],
'Kfold-Mean-Accuracy':[results.mean()],
'Kfold-StdDeviation':[results.std()]}))
In [583]:
modelComp
Out[583]:
In [ ]:
Attempting Random Forest with PCA to see the impact¶
In [501]:
X=df_scaled.iloc[:,0:8]
y = df_scaled.iloc[:,8]
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.3, random_state = seed)
from sklearn.decomposition import PCA
pca = PCA(8)# Initialize PCA object
pca.fit(X_train)
Out[501]:
In [502]:
pca.explained_variance_
Out[502]:
In [503]:
pca = PCA(n_components=6)
pca.fit(X_train)
X_train_pca = pca.transform(X_train) # PCs for the train data
X_test_pca = pca.transform(X_test) # PCs for the test data
X_train_pca.shape, X_test_pca.shape
Out[503]:
In [504]:
rf = RandomForestRegressor(bootstrap = True,
max_depth= 15,
max_features= 'sqrt',
min_samples_leaf= 1,
n_estimators= 100)
rf.fit(X_train_pca, y_train)
y_pred = rf.predict(X_test_pca)
# Train data accuracy
print('Performance on training data using DT:',rf.score(X_train_pca,y_train))
# test data accuracy
print('Performance on testing data using DT:',rf.score(X_test_pca,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
model1 = rf
results = cross_val_score(model1, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
- With PCA implemetation and our existing data set (with selective attributes), we can see similar impact
- we will proceed with existing data set and not with PCA implemented data set for better understanding of relations
Gradient Boosting Regressor¶
In [594]:
X=df_scaled.iloc[:,0:8]
y = df_scaled.iloc[:,8]
#Removing less contributing features
X=X.drop(['ash','coarseagg','fineagg'],axis=1)
# Split X and y into training and test set in 70:30 ratio
X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.3, random_state = seed)
In [506]:
model=GradientBoostingRegressor()
model.fit(X_train, y_train)
Out[506]:
In [507]:
y_pred = model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
model1 = model
results = cross_val_score(model1, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
- With Gradiant boost regressor we can notice the train and test set accuracies are close by compared to other models with default setting itself
- We will try to get better close by numbers for train and Test to make it more generalized with out loosing much of test accuracies.
- Cross val score is also better at accuracy of 85%
we can take the following approach:¶
- Choose a relatively high learning rate. Generally the default value of 0.1 works but somewhere between 0.05 to 0.2 should work for different problems
- Determine the optimum number of trees for this learning rate. This should range around 5-10 as we noticed earlier dring decission tree as well. Remember to choose a value on which your system can work fairly fast. This is because it will be used for testing various scenarios and determining the tree parameters.
- Tune tree-specific parameters for decided learning rate and number of trees. Note that we can choose different parameters to define a tree .
- Lower the learning rate and increase the estimators proportionally to get more robust models.
In [ ]:
In [508]:
parameters = {
"loss":['ls', 'lad', 'huber', 'quantile'],
"learning_rate": [0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2],
#"min_samples_split": np.linspace(0.1, 0.5, 5),
#"min_samples_leaf": np.linspace(0.1, 0.5, 5),
"max_depth":[3,5,8],
"max_features":["log2","sqrt"],
"criterion": ["friedman_mse", "mae"],
"subsample":[0.5, 0.618, 0.8, 0.85, 0.9, 0.95, 1.0],
"n_estimators":[10]
}
clf = GridSearchCV(estimator = GradientBoostingRegressor(),
param_grid = parameters,
cv = 5,
verbose = 2,
n_jobs= 4)
clf.fit(X_train, y_train)
clf.best_params_
Out[508]:
In [510]:
#Splitted the param into two reduce the time
parameters = {
"loss":['ls'],
"learning_rate": [0.2],
"min_samples_split": np.linspace(0.1, 0.5, 5),
"min_samples_leaf": np.linspace(0.1, 0.5, 5),
"max_depth":[8],
"max_features":["log2"],
"criterion": ["friedman_mse"],
"subsample":[0.9],
"n_estimators":[10,100]
}
clf = GridSearchCV(estimator = GradientBoostingRegressor(),
param_grid = parameters,
cv = 5,
verbose = 2,
n_jobs= 4)
clf.fit(X_train, y_train)
clf.best_params_
Out[510]:
model=GradientBoostingRegressor(criterion= 'mae', learning_rate= 0.2, loss= 'huber', max_depth= 5, max_features= 'sqrt', n_estimators= 100, subsample= 0.9, min_samples_leaf=0.1, min_samples_split= 0.2,) model.fit(X_train, y_train)
In [584]:
model=GradientBoostingRegressor(criterion= 'friedman_mse',
learning_rate= 0.2,
loss= 'ls',
max_depth= 8,
max_features= 'log2',
n_estimators= 100,
subsample= 0.9,
min_samples_leaf=0.1,
min_samples_split= 0.2,)
model.fit(X_train, y_train)
Out[584]:
In [585]:
y_pred = model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
#model1 = model
results = cross_val_score(model, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
modelComp=modelComp.append(pd.DataFrame({'Model':['GradientBoostingRegressor'],
'Train Accuracy':[model.score(X_train,y_train)],
'Test Accuracy':[model.score(X_test , y_test)],
'Kfold-Mean-Accuracy':[results.mean()],
'Kfold-StdDeviation':[results.std()]}))
- Traing and test data sets accuracies are very close around 90% and the cross val score is as well close to 85%
- This model is more geenralised as the train and test accuracies are very close at 92% and 89% respectively
In [586]:
modelComp
Out[586]:
In [ ]:
AdaBoostRegressor¶
In [597]:
model=AdaBoostRegressor()
model.fit(X_train, y_train)
Out[597]:
In [598]:
y_pred = model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
model1 = model
results = cross_val_score(model1, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
In [ ]:
#Custom prameters
In [603]:
model=AdaBoostRegressor(n_estimators=100,
learning_rate=1,
loss='linear')
model.fit(X_train, y_train)
Out[603]:
In [604]:
y_pred = model.predict(X_test)
# Train data accuracy
print('Performance on training data using DT:',model.score(X_train,y_train))
# test data accuracy
print('Performance on testing data using DT:',model.score(X_test,y_test))
#r2 score
acc_DT=metrics.r2_score(y_test, y_pred)
print('Accuracy Test: ',acc_DT)
print('MSE: ',metrics.mean_squared_error(y_test, y_pred))
kfold = KFold(n_splits=num_folds, random_state=seed, shuffle=True)
#kfold = KFold(n_splits=num_folds)
model1 = model
results = cross_val_score(model1, X_train, y_train, cv=kfold)
print(results)
print("\n Average model Accuracy: %.3f%% with std. dev - (%.3f%%)" % (results.mean()*100.0, results.std()*100.0))
- AdaBoost Regressor Traing and Test set accuracies are around 80%, so the model is a not a over fit or under fit model.
- But the cross val score is low compared to other regreesors
- We were able to sqeeze out little extra gain in both train and test accuracies but its still not comparable to other regressors.
In [605]:
modelComp
Out[605]:
- Comparing these selective details, we can notice the train and test accuracies are closest for gradient boost followed by Random forest regressor
- Keeping in mind about overfitting modles we are tilting towrds GBM it has showed similar test accuracy and train accuracies were also close by.
- So far we can consider Gradient boost as its been more generalized and its Train/Test/cross-val accuracies are nearby compared to other models
Bootstrap Sampling and Confidence Interval¶
In [395]:
from sklearn.utils import resample
from sklearn.metrics import accuracy_score
from matplotlib import pyplot
data = X.join(y)
values = data.values
In [ ]:
In [401]:
#Random Forest
# Number of bootstrap samples to create
n_iterations = 1000
# size of a bootstrap sample
n_size = int(len(data) * 0.5)
# run bootstrap
# empty list that will hold the scores for each bootstrap iteration
stats = list()
for i in range(n_iterations):
# prepare train and test sets
train = resample(values, n_samples=n_size) # Sampling with replacement
test = np.array([x for x in values if x.tolist() not in train.tolist()]) # picking rest of the data not considered in sample
# fit model
model = RandomForestRegressor(
bootstrap = True,
max_depth= 8,
max_features= 'sqrt',
min_samples_leaf= 1,
n_estimators= 100)
# fit against independent variables and corresponding target values
model.fit(train[:,:-1], train[:,-1])
# Take the target column for all rows in test set
y_test = test[:,-1]
# evaluate model
# predict based on independent variables in the test data
score = model.score(test[:, :-1] , y_test)
predictions = model.predict(test[:, :-1])
stats.append(score)
In [402]:
# plot scores
pyplot.hist(stats)
pyplot.show()
# confidence intervals
alpha = 0.95 # for 95% confidence
p = ((1.0-alpha)/2.0) * 100 # tail regions on right and left .25 on each side indicated by P value (border)
lower = max(0.0, np.percentile(stats, p))
p = (alpha+((1.0-alpha)/2.0)) * 100
upper = min(1.0, np.percentile(stats, p))
print('%.1f confidence interval %.1f%% and %.1f%%' % (alpha*100, lower*100, upper*100))
In [ ]:
In [398]:
#Gradient Boost
# Number of bootstrap samples to create
n_iterations = 1000
# size of a bootstrap sample
n_size = int(len(data) * 0.5)
# run bootstrap
# empty list that will hold the scores for each bootstrap iteration
stats = list()
for i in range(n_iterations):
# prepare train and test sets
train = resample(values, n_samples=n_size) # Sampling with replacement
test = np.array([x for x in values if x.tolist() not in train.tolist()]) # picking rest of the data not considered in sample
# fit model
model = GradientBoostingRegressor(criterion= 'mae',
learning_rate= 0.2,
loss= 'huber',
max_depth= 5,
max_features= 'sqrt',
n_estimators= 100,
subsample= 0.9,
min_samples_leaf=0.1,
min_samples_split= 0.2,)
# fit against independent variables and corresponding target values
model.fit(train[:,:-1], train[:,-1])
# Take the target column for all rows in test set
y_test = test[:,-1]
# evaluate model
# predict based on independent variables in the test data
score = model.score(test[:, :-1] , y_test)
predictions = model.predict(test[:, :-1])
stats.append(score)
In [400]:
# plot scores
pyplot.hist(stats)
pyplot.show()
# confidence intervals
alpha = 0.95 # for 95% confidence
p = ((1.0-alpha)/2.0) * 100 # tail regions on right and left .25 on each side indicated by P value (border)
lower = max(0.0, np.percentile(stats, p))
p = (alpha+((1.0-alpha)/2.0)) * 100
upper = min(1.0, np.percentile(stats, p))
print('%.1f confidence interval %.1f%% and %.1f%%' % (alpha*100, lower*100, upper*100))
In [ ]:
- For Gradient Boost - 95.0 confidence interval 82.6% and 87.5%
- For Random Forest - 95.0 confidence interval 81.6% and 86.6%
- there is s slight advantage for Gradient Boost model compared to Random Forest
- Comparing the model performance range at 95% confidence we can notice its a tough decission to decide between the models solely comparing the number. We have to weigh in the model design and what need to be considered for the data and complexity.
- RF basically has only one hyperparameter to set: the number of features to randomly select at each node. However there is a rule-of-thumb to use the square root of the number of total features which works pretty well in most cases but need to look upon case by case. On the other hand, GBMs have several hyperparameters that include the number of trees, the depth (or number of leaves), and the shrinkage (or learning rate).
- There is one fundamental difference in performance between the two that may force you to choose Random Forests over Gradient Boosted Machines (GBMs). That is, Random Forests can be easily deployed in a distributed fashion due to the fact that they can run in parallel, whereas Gradient Boosted Machines only run trial after trial.
- Considering our scenario we have limited set of data hence the the cons for parallelism is not affecting us heavily. Along with that our data has very less linear relations and GBM approach of multiple iteration for adding weight might give us added advantage over RF.
Reasons for choosing GBM¶
- RF overfit a sample of the training data and then reduces the overfit by simple averaging the predictors. But GBM repeatedly train trees or the residuals of the previous predictors.
- RF is easy to use (less tuning parameters). we can blindly apply RF and can get decent performance with a little chance of overfit but without cross validation GBM is useless. GBM need much care to setup. As in GM we can tune the hyperparameters like no of trees, depth, learning rate so the prediction and performance is better than the Random forest.
- With continuos learning and exposure to data we can fine tune GBM more accurately and in generalized form
In [ ]: